#################################################################################################
#################################################################################################
####Appendix: Non-state Atrocities in Capital Cities - ZINB Nighttime Light Robustness Models####
#################################################################################################
#################################################################################################
library(MASS) 
library(pscl)
library(foreign)
library(stargazer)
library(mvtnorm)
library(plyr)
library(glmmADMB)
library(MCMCglmm)

# Set working library
setwd("~/Data/Global Analysis/")
# Read in main data
main.data <- read.dta("full.grid.upd.dta")

############################
###Global sample analysis###
############################
###Model A-XV (Baseline model)
at.zinb.ai.xv <- zeroinfl(incidentnonstatefull ~ lag_Capital + laglnNL_sum + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                      + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                        year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 
                      |laglnNL_sum + lagurban + loglagttime + loglagcellarea, 
                      dist = "negbin", data = main.data)
summary(at.zinb.ai.xv)
AIC(at.zinb.ai.xv)

###Model A-XVI (Full model)
at.zinb.ai.xvi <- zeroinfl(incidentnonstatefull ~ lag_Capital + laglnNL_sum +lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                      + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                        year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 
                      |laglnNL_sum + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                      + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                      dist = "negbin", data = main.data)
summary(at.zinb.ai.xvi)
AIC(at.zinb.ai.xvi)

#################################
###Urban cells sample analysis###
#################################
#Subset data
u.main.data <- main.data[ which(main.data$lagurban > 0),]
###Model A-XVII (Baseline model)
at.zinb.ai.xvii <- zeroinfl(incidentnonstatefull ~ lag_Capital + laglnNL_sum + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                          year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 
                        |laglnNL_sum + lagurban + loglagttime + loglagcellarea, 
                        dist = "negbin", data = u.main.data)
summary(at.zinb.ai.xvii)
AIC(at.zinb.ai.xvii)

###Model A-XVIII (Full model)
at.zinb.ai.xviii <- zeroinfl(incidentnonstatefull ~ lag_Capital + laglnNL_sum + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                          year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 
                        |laglnNL_sum + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                        + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                        dist = "negbin", data = u.main.data)
summary(at.zinb.ai.xviii)
AIC(at.zinb.ai.xviii)

#####Export to LaTex
stargazer(at.zinb.ai.xv, at.zinb.ai.xvi, at.zinb.ai.xvii, at.zinb.ai.xviii)
stargazer(at.zinb.ai.xv, at.zinb.ai.xvi, at.zinb.ai.xvii, at.zinb.ai.xviii, zero.component = TRUE)
AIC(at.zinb.ai.xv, at.zinb.ai.xvi, at.zinb.ai.xvii, at.zinb.ai.xviii)
